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ABSTRACT 

Context. In spite of enormous progress and brilliant achievements in cosmological simulations, they still lack numerical 
resolution or physical processes to simulate dwarf galaxies in sufficient details. Accurate numerical simulations of indi¬ 
vidual dwarf galaxies are thus still in demand. 

Aims. We aim at improving available numerical techniques to simulate individual dwarf galaxies. In particular, we aim 
at (i) studying in detail the coupling between stars and gas in a galaxy, exploiting the so-called stellar hydrodynamical 
approach, and {ii) studying for the first time the chemo-dynamical evolution of individual galaxies starting from self- 
consistently calculated initial gas distributions. 

Methods. We present a novel chemo-dynamical code for studying the evolution of individual dwarf galaxies. In this code, 
the dynamics of gas is computed using the usual hydrodynamics equations, while the dynamics of stars is described by 
the stellar hydrodynamics approach, which solves for the first three moments of the collisionless Boltzmann equation. 
The feedback from stellar winds and dying stars is followed in detail. In particular, a novel and detailed approach has 
been developed to trace the aging of various stellar populations, which enables an accurate calculation of the stellar 
feedback depending on the stellar age. The code has been accurately benchmarked, allowing us to provide a recipe for 
improving the code performance on the Sedov test problem. 

Results. We build initial equilibrium models of dwarf galaxies that take gas self-gravity into account and present dif¬ 
ferent levels of rotational support. Models with high rotational support (and hence high degrees of flattening) develop 
prominent bipolar outflows; a newly-born stellar population in these models is preferentially concentrated to the galactic 
midplane. Models with little rotational support blow away a large fraction of the gas and the resulting stellar distri¬ 
bution is extended and diffuse. Models that start from non-self-gravitating initial equilibrium configurations, evolve at 
a much slower pace owing to lower initial gas density and hence lower star formation rate. The stellar dynamics turns 
out to be a crucial aspect of galaxy evolution. If we artificially suppress stellar dynamics, supernova explosions occur in 
a medium heated and diluted by the previous activity of stellar winds, thus artificially enhancing the stellar feedback. 
Conclusions. The stellar hydrodynamics approach presents a promising tool for studying numerically the coupled evo¬ 
lution of gas and stars in dwarf galaxies. 

Key words. Galaxies: abundances - Galaxies: dwarf - Galaxies: evolution - Galaxies: ISM - Galaxies: kinematics and 
dynamics 


1. Introduction 

Cosmological simulations are reaching nowadays extraordi¬ 
nary levels of complexity and sophistication. Based on the 
so-called standard model of cosmology (SMoC), these state- 
of-the-art simulations are able to repro duce the ma i n fea- 


tures of observed large galaxies (see e.g 

iBrook et ah. 2012 

Zem 

p et ah. 20121 

iHoDkins et alJ.l201,‘i 

Voffelsberser et al. 

2014 

: Ferret et al 

.1. 120l4 ISchave et al.l 

2014ll. althousfh se- 


vere disagreement s with observations still persist (see e.g. 
iGenel et al.l . 1201^ . However, the resolution of these simu¬ 
lations still does not allow to treat small-scale physics in 
dwarf galaxies (DCs). 

Unfortunately, most of the weaknesses and prob¬ 
lems shown by the SMoC manifest the mselves i n the 
range of masses typi c al for DCs ( see e.g. iKroupal 120121: 
iFamaev fc McGanghL 120121 l201.ll for recent reviews). 
Besides full-blown cosmological simulations, it is thus 
necessary to run simulations of individual galaxies or 
small groups of galaxies. This enables studying physi¬ 


cal processes in more numerical resolution and address 
the problems of the SMoC. The simulation of individ¬ 
ual galaxies is therefore still extremely relevant. Also in 
this field very remarkable level of accuracy have been 
reached (see e. g. iH onkins et al.l. I ' 

Roskar et al.l . 1201,11 iMinchev et 


MM 

: all l2 


■ .Renaud et al 

201 . 11 : 


2013 


- . — ■ _ _ _ — ■ _ iTevssier et al._ 

2013L among others). The simulation of individual galax¬ 


ies can also be used to parameterise feedback effects that 
can not be treated i n detai l in c osmological simulations 
(jCreasev et al.l . l2013l |R ecc hil. [2nil . 


An obvious drawback of simulating individual galaxies 
is that it is not clear what initial conditions one should 
adopt. A common strategy is to consider a rotating, isother¬ 
mal gas in equilibrium with the potential generated by a 


fixed distribution of stars and/or of dark matter (DM), 
but not with the potential generated by the gas itself 

(iTomisaka & Ikenchi 

. 19881: ISilich & Tenorio-Taglel Il998l 

iMac-Low & Ferraral 

1999HRecchi & Hensler. 20131). A tvp- 

ical justihcation for 

neglecting gas self-gravity in DC sim- 
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ulations is that the mass budget of these objects should 
be dominated by dar k matter. In our previous paper 
(|Vorobvov et al.l . l2012f) . we have explained in detail why 
this approach may produce unrealistic initial gas distribu¬ 
tions and considered also the gas self-gravity when building 
equilibrium initial conditions. 

As known from st ellar structure studies (see e.g. 
lOstriker fc Mark! . Il968h . the inclusion of gas self-gravity 
leads to implicit equations describing the initial equilib¬ 
rium conhguration, which must be solved iteratively. We 
have done that for different halo masses, different degrees 
of rotational support and different initial temperatures and 
we have calculated the amount of gas prone to star forma¬ 
tion for each of these models. Our star formation criteria 
were based on the Toomre instability criterium and on the 
Kennicutt-Schmidt empirical correlation between gas sur¬ 
face density and sta r form ation rate (SFR) per unit area 
(see IVorobvov et al.l . l2012l for more details). For models 
that satisfy both star formation criteria, the underlying as¬ 
sumption was that the initial configuration is marginally 
stable to star formation and an external perturbing agent 
can trigger star formation in the corresponding parts of the 
galaxy. 


2. Simulation methodology 

Our model DCs consist of a gas disk, stellar component, 
and fixed dark matter halo. The gas disk is a mixture of 
nine chemical elements (H, He, C, O, N, Ne, Mg, Si, and 
Fe), with the initial abundance of heavy elements set by the 
adopted initial metallicity. The stellar component consists 
of new stars born in the course of numerical simulations 
and has no pre-existing component. The DM halo has a 
spherically symmetric form and contributes only to the to¬ 
tal gravitational potential of the system. 

We describe the time evolution of our model galaxies 
using the coupled system of gas and stellar hydrodynam¬ 
ics equations complemented with phase transformations be¬ 
tween stars and gas, chemical enrichment by supernova ex¬ 
plosions, low- and intermediate-mass stars and also star 
formation feedback. We assume that individual chemical 
elements are well coupled dynamically with the bulk mo¬ 
tion of the gas, which means that we need to solve only 
for the continuity equation for every chemical element with 
mass density pi and need not to solve for the dynamics of 
every element separately. Below, we provide a brief expla¬ 
nation for the key equations used to evolve our system in 
time. 


The subsequent logical step is to use these marginally 
stable configurations as initial conditions for full-blown hy- 
drodynamical simulations. For this purpose, we have de¬ 
veloped a detailed 2+1 dimensional hydrodynamical code 
in cylindrical coordinates with assumed axial symmetry, 
which treats both the stellar and gaseous components as 
well as the phase transitions between them. Of course, a 
real galaxy may depart from such axial symmetry. However, 
since our initial conditions are axially symmetric and we do 
not explicitly consider environmental effects in this paper, 
the departures from axial symmetry are not expected to be 
significant. 


This is the first of a series of papers devoted to numer¬ 
ically studying DGs with self-consistent initial conditions. 
Here, we explain in detail the basic strategy and the em¬ 
ployed numerical hydrodynamical methods, which we ex¬ 
tensively benchmark. Furthermore, we run a few represen¬ 
tative models to show the overall early evolution of a typ¬ 
ical DG. We pay particular attention to the co-evolution 
of stars and gas in our simulations. We follow the mo¬ 
tion of stars by means of the so-called stellar hydrody¬ 
namical approach (iTheis et al.l. 1992^ Samland et al.L 1997; 


Vorobyov fc Theil l2006t iMitchell et al.l . 120131 : iKnlikov . 

201411 . according to which the stellar component of a galaxy 
can be treated as a fluid. We develop this approach further 
and describe how stars and gas can be coupled by means 
of star formation and stellar feedback processes. 


The paper is organised as follows. In section [5] we 
describe the simulation methodology, paying specific at¬ 
tention to coupling the stars and gas via star formation 
and stellar feedback. Section |3] contains the gas cooling 
and heating rates employed in our modeling. The solution 
method of our stellar hydrodynamics equations is described 
in Section 0] The evolution of several representative mod¬ 
els of DGs is studied in Section 0 Model caveats are dis¬ 
cussed in Section [6] The main conclusions are summarized 
in Section [7] and we benchmark our code against a suite of 
test problems in the Appendix. 


2.1. Gas hydrodynamics 

The dynamics of gas is modelled using the usual hydro¬ 
dynamics equations for the gas density of pg, momentum 
PgVg, and internal energy density e complemented with the 
rates of mass, momentum and energy transfer between the 
gas and stellar components and also with the continuity 
equations for the mass density of every chemical element p; 

^ + V-ipgVg)=S-V, (1) 

^ + V • ip,vg) =S^- V\ (2) 

d 

^ (Pg^g) + V • {pgVg ® Vg) = -VP + Pg£/g, + 

+ v^S - VgV, (3) 

^ + V-(e'i;g) = -PV-t;g-A + r + r,-e—. (4) 

ot Pg 

Here, Vg and Vs are the gas and stellar velocities, P = 
(7 — l)e is the gas pressure linked to the internal energy 
density via the ideal equation of state with the ratio of 
specific heats 7 = 5/3 and g is the gravitational accel¬ 
eration due to gas, stars, and DM halo. The gas cooling 
A due to line and continuum emission and gas heating F 
due to cosmic rays and small dust grains are explained in 
detail in Section [3] The gas heating due to stellar feedback 
F* is calculated as the sum of the contributions from SNH 
and SNIa explosions (Fsnii and FsNia, respectively), and 
also from the stellar winds Fsw The last term on the r.h.s. 
of equation o accounts for the decrease in the internal 
energy due to formation of stars from the gas phase. We 
note that the quantity PgVg (g) Vg is the symmetric dyadic 
(a rank-two tensor) calculated according to the usual rules 
and V • {pgVg 0 Vg) is the divergence of a rank-two tensor. 

The source term S denotes the phase transformation 
of stars into the gas phase and is defined as the sum of 
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the mass release rates (per unit volume) of all individual 
elements, S = The mass release rate of a particular 

element i is defined as 

*5* = ‘^SNII + ‘^SNIa + '5*) (5) 

where and 5® are the contributions due to 

SNII, SNIa, and also due to low- and intermediate mass 
stars, respectively. The sink term 1) describes the phase 
transformation of gas into stars according to the adopted 
star formation law. The detailed explanation of all relevant 
rates is provided in Section [2T1 


2.2. Stellar hydrodynamics 

The time evolution of the stellar component is computed 
using the Boltzmann moment equ ation approach intro¬ 
duced by iBurkert fc Henslerl (Il988fland further devel oped 
in application to stellar d isks in ISamland et al.l (|l997fl and 
IVorobvov fc Theii (l2006ll . In this approach, the dynamics 
of stars is described by the first three moments of the col¬ 
lisionless Boltzmann moment equations 

^ + V ■{p,v,)=V-S, (6) 

d 

— (psVs) + V • 0 PsVs) = -V • n 4 - Ps5gr - 

- VsS + VgV, (7) 


-f V • (n^-Vs) = -liik : {Vvs)jk - Ujk ■■ {Vvs)ik - 

- alS + alV, (8) 

where ps = rn J fd^u is the stellar mass density, Vs = 
Ps^rn f fud^u is the mean stellar velocity, and 11 is the 
stress tensor with six components 11^ = PsCrfj. The stel¬ 
lar velocity dispersions are defined as crfj = p~^m f f (ui — 
Vi)(uj — Vj)d^u. The function / is the distribution func¬ 
tion of stars in the six dimensional position-velocity phase 
space / = f{t,x,u) and m is the average mass of a star. 
The quantity is the gradient of a ve ctor, the covariant 
expres sion for which is provided, e.g., in IStone fc NormanI 
(1199211 and : (Vt)s)jfe is the convolution (over index k) 
of two rank-two tensors. The quantity ct* represents a typi¬ 
cal stellar velocity dispersion of a new-born cluster of stars, 
for which we take a fiducial value of 5.0 km s“^. 


2.3. Massive stars 

Massive stars that end up their life with SNII explosions 
are main contributors to the energy budget and chemical 
composition of the interstellar medium. Therefore, it is im¬ 
portant to follow their evolution as accurate as possible. We 
improve the description of the stellar component by solv¬ 
ing for a separate continuity equation for the mass volume 
density Pg of stars with mass >8.0 Mq 

^ + V•(p^ns)=/hI?-5^ (9) 

where /h is the fraction of stars with mass > 8.0 Mq calcu¬ 
lated according to the adopted IMF (see equation [T3]) and 


5h is the rate of death of massive stars. The massive star 
subcomponent has the same dynamical properties as the 
rest of the stars, an assumption that we plan to relax in 
the future. 


2.4. Intermediate-mass stars 


Stars in the (1.0 — 8.0) Mq mass range produce the ma¬ 
jority of carbon and nitrogen. Moreover, binary stars in 
this mass range can explode as Type la supernovae and are 
thus important sources of iron. Therefore, we also solve for 
a separate continuity equation for the mass volume density 
p™ of intermediate-mass stars 



+ V • (p“^.) 


fmT) - 5” 


( 10 ) 


where /m is the fraction of stars with mass 1.0 Mq < m < 
8.0 Mq calculated according to the adopted IMF and 5i„ 
is the rate of death of intermediate-mass stars. We note 
that the stellar evolution of low-mass stars can be neglected 
during the galactic evolution. They only contribute to the 
gravitational potential where their density can be deduced 
from Ps (the total stellar density), p^, and p™. 


2.5. Gravity of gas, stars, and DM halo 

The gravitational potential of gas and stars is calculated 
self-consistently by solving for the Poisson equation 

V^<i)g+s = 47rG (pg-I-Ps) ■ (11) 

To accelerate the simulations, we solve for Equation (HU 
only when the relative error in the currently stored gravi¬ 
tational potential com pared to the current den sity distribu¬ 
tion exceeds 10“® ('see lStone &: NormanI . 1199^ for details). 

The gravitational acceleration due to a spherically 
symmetric DM halo is calculated as explained in 
IVorobvov et al.l (l2012f) . Two options are available in the 
code: a cored isothermal DM halo profile and a Navarro- 
Frenk-White (or cuspy) profile. In this paper, we use the 
cored isothermal DM profile. The total gravitational accel¬ 
eration is calculated as = Oh ~ ^4>g+s. 


2.6. Mean stellar age 


In order to calculate the mass and energy release rates by 
supernovae and dying low-mass stars one needs to know 
the age of stars. In the stellar hydrodynamics approach, 
stars are a mixture of populations with different ages. 
Nevertheless, one may define the mean age of stars t*, 
which should obey the following rules: 


1. The bulk motion of stars, such as compression or rar¬ 
efaction, should not affect the value of r*. This can 
be achieved by solving for the Lagrangian or comoving 
equation for r* 


cff* 


= 0 


( 12 ) 


2. A newly-born stellar population should rejuvenate the 
pre-existing one. This is achieved by updating the mean 
age every timestep in every grid cell using the following 
equation 


MgT* -I- AMs At 

Ms -f AMs 


(13) 
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Radial distance Radial distance 


Fig. 1. Star volume density (left) and mean stellar age 
(right) as a function of time during the gravitational col¬ 
lapse of a stellar sphere of unit radius. The elapsed time is 
indicated in the left panel. Star formation is turned off for 
this test problem. The dash-dot-doted lines shows the mean 
stellar age when calculated using the Eulerian form to ac¬ 
count for bulk motions of stars instead of the Lagrangian 
form (see text for more detail.) 


where Mg is the stellar mass in a specific cell and AMg 
the mass of stars born in the same cell during timestep 
At. Equation (IT^ satisfies the expected asymptotic be¬ 
haviour of the stellar age, which becomes small (r^-^'' —?► 
At) during a massive instantaneous burst (AMg —oo), 
but remains almost unchanged (tJ'' —?► r*) in the qui¬ 
escent phase AMg —>■ 0. 

3. The mean stellar age should uniformly increase with 
time, reflecting the overall aging of the system. This 
is done by adding the current timestep At to the mean 
stellar age in every grid cell. Thus, in the absence of star 
formation and pre-existing old stellar population, the 
mean stellar age is equivalent to the current evolution 
time t. 

In practice, we first solve for Equation (HU to account 
for bulk motions, then we update the mean stellar age in 
every grid cell according to Equation (fT51) to account for 
star formation, and finally we add the current timestep At 
to the mean stellar age in every grid cell to take into ac¬ 
count the aging of stellar populations. Since we follow the 
evolution of both the massive and intermediate-mass star 
subcomponents, their mean ages and t™, respectively, 
need to be calculated using the same procedure. 

Figure [T] shows the time evolution of the stellar den¬ 
sity ps (left panel) and mean stellar age r* (right panel) 
in an idealized test problem involving the gravitational col¬ 
lapse of a stellar sphere with unit radius. Initially, ps and 
r* are set to unity inside the sphere and to zero outside the 
sphere (dash-dotted lines). The stellar velocity dispersion 
is negligible everywhere and the star formation rate is set 
to zero. As the collapse proceeds, stellar density shows the 
expected behaviour with a central plateau shrinking in size 
and growing in density. At the same time the mean stel¬ 
lar age increases in concordance with time, as expected in 
the absence of star formation. In particular, the difference 
between the values of r* inside and outside the sphere is al¬ 
ways equal to unity, as was set initially. We emphasize here 
the need for using the Lagrangian equation for r*. If the 
Eulerian form -|- V • {vsf^) = 0 is used instead, the 

above test fails. Indeed, the dash-dot-dotted line shows the 
mean stellar age at t = 0.4 calculated using the Eulerian 
form. Evidently, the difference between the values of r* in¬ 


side and outside the sphere is now much greater than unity, 
indicating a spurious aging of the stellar populations inside 
the sphercQ. 

Finally, we note that the concept of the mean stellar 
age has its limitations. For instance, in the case of a con¬ 
stant SFR and steady galaxy conhguration, approaches 
a constant value after a few tens of Myr, meaning that the 
supernova rates and mass return rates (see equations (ITTl) 
and HU)) will be exclusively determined by stars whose life¬ 
time is equal to r^. This bias towards stars with single age 
(and mass) diminishes, when a time-varying SFR is present 
or stellar motions are taken into account. Nevertheless, a 
more sophisticated approach that takes into account a pos¬ 
sible age spread around the mean value is desirable. 


2.7. Source and sink terms 


In order to calculate the rates at which stars return their 
mass (including newly synthesized elements) into the gas 
phase, we need to make assumptions about the initial mass 
function (IMF), stellar lifetimes and nucleosyn thesis rates. 

We adopt the Kroupa IMF (iKroupal l200Ill of the form 


<(>(M,) 


AM-1'3 ifM*<O.5M0 
B M-2-3 if M* > 0.5 Mq , 


(14) 


where the normalization constants A and B are calculated 
as described in Sections [2.7.11 and 12.7.21 

The stellar l ifetim es are taken from 

iPadovani fc Matteuccil (|1993ll 


t(M*) 


1.2 M-i ®5 + 0.003 Gyr if M* > 7.45 M© , . 
10/(M.) Gyr if M* < 7.45 Mq , 


where 

f{M.) 


|o.334- v'1.79 - 0.2232 x (7.764 - log(M*))| 
0.1116. 


(16) 


In equation (I15L we have adjusted the transitional value 
for the stellar mass (7.45 M©) in order to obtain a smooth 
time derivative dM^/dr. 


2.7.1. Supernova type II rates 

We assume that all stars in the [8.0 M© : mass 

range end their life cycle as type II supernovae (SNell), 
where is the upper cut-off mass of the IMF set to 

100 Mq in our model. The SNII rate, i.e., the number of 
SNell per unit time, is then defined as 


l^SNII = 


ifrje [t(100M©):t(8.0M©)], 
ifr (0 [r(lOOM0):T(8.OM©)], 


(17) 


0 


where is the mass of a massive star that is about 

to die as SNII calculated by inverting equation (ITH1) . This 
is the SNII rate following an instantaneous burst of star 
formation. We therefore assume here that each star forma¬ 
tion event generates a single stellar population (SSP). The 

^ It is still possible to use the Eulerian form for the 
product psV* instead of just r*. We leave the exact 
derivation of the pertinent equations for a followup 
study. 
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superposition of SNII rates coming from many SSPs will 
then approximate the SNII rate in the presence of a vari¬ 
able SFR. The derivative is calculated by invert¬ 

ing equation (HH) and then differentiating it with respect 
to time. 

The IMF in equation m has to be normalized accord¬ 
ing to the total mass of massive stars in each com¬ 

putational cell 

J (18) 

8.0 Mq 

We note that the upper limit in the integral is not fixed but 
rather depends on the mean stellar age of massive stars 
in a given computational cell. This is done to take into 
account the fact that massive stars with lifetimes smaller 
than T^, and hence with mass greater than must 

have exploded as SNell. Excluding those stars ensures a 
correct calculation of the stellar mass spectrum at any time 
instanlQ. Since and are time-varying quantities, 

the normalization constants in the IMF need to be updated 
every time step. 

The mass return rate (per unit volume) by SNII explo¬ 
sions can now be defined as 

'^SNII = ‘^SNII PSNIl[-^*(''■*)]) (19) 

where Psnii[-^*( 1 ^*)] mass (per unit volume) released 

by a SNII of mass M* in the form of a specific element 
i. The values of p® for nine considered elements (H, He, 
C, O, N, Ne, Mg, Si, and Fe) and four stellar metallicities 
(Z^ = 10 ~^, 10 ~^, 10 ~^, an d 1.0 times solar) are taken from 
IWooslev fc WeavCT () 199511 . For stellar masses greater than 
40 Mq, the yields are set constant to those of a 40 Mq star. 

Finally, we assume that each supernova releases esN x 
10 ®^ ergs in the form of thermal energy and define the rate 
of energy density released by SNell as 

Fsnii = lO^^esNl^SNii V~^, (20) 

where V is the injection volume specified by the actual size 
of our numerical grid and esN is the ejection efficiency set 
to 1.0 in this work (see Section [5] for discussion.). 

2.7.2. Supernova type la rates 

We derive the rates of Type la supernovae follow i ng th e 
ideas and notations laid out in iMatteucci fc Recchil (l2001ll . 
The number of binary systems (per unit mass of the binary 
system Mb) that are capable of producing a SNIa explosion 
can be calculated as 


defined as 

Mb,i = max[ 2 M 2 (T“),MB.mi„] ( 22 ) 

Mb,2 = -MB,max + M2 (t™). ( 23 ) 

The maximum and minimum masses of a binary system 
that is capable of producing a SNIa explosion are set to 
Mb, min = 3.0 Mq and Mb, max = 16.0 Mq, respectively, 
and the maximum mass of the secondary is M 2 = 8.0 Mq. 
The quantity /(M 2 (t“)/Mb) in equation (1^ is the dis¬ 
tribution function of the M 2 (t™)/Mb ratio and is defined 
as 



and 7 is set to 0.5. For this value of 7 , the normalization 
constant that bettter reproduces the observed [a / Fe] ra - 
tios in dwarf galaxies is = 0.06 (iRecchi et al.l . l2009ll . 
therefore we adopt this value in the present work. 

By analogy to the massive star subcomponent, the IMF 
in equation m has to be normalized according to the to¬ 
tal mass of intermediate-mass stars Mg™^^ in each compu¬ 
tational cell 

(t”) 

J (j){M,)M^ dM,, = Mg“of (25) 

1.0 Mq 

The SNIa rate, i.e., the number of SNIa explosions 
per unit time, can now be calculated by analogy to equa¬ 
tion (fT71) as 

f iVsNIa[M2(r“))]^^^^^^ 

"^SNia = < if G [r(8.0 Mq) : r(M 2 ,min)], (26) 

[0 if t“ 9 [r(8.0 Mq) : r(M 2 ,min)], 

where M 2 , min =1.0 Mq is the minimum mass of the sec¬ 
ondary. 

The mass return rate (per unit volume) by SNIa explo¬ 
sions can now be expressed as 

‘^SNIa = ’^SNIa PsNIa[-^*(''’™)]i (27) 

where PsNia[-^*(l^D] i® Ih® mass (per unit volume) released 
by a SNIa of mass M* in the form of a specific element i. 
The values of for six elements (C , O, N, Mg, Si, and Fe) 
are taken from lNomoto et H] (Il984f) . We neglect the metal- 
licity dependence of the SNIa yields. The energy release rate 
by SNIa (FsNia) is defined in analogy to Equation (l20l) . with 
T^snii being substituted by TT-sNia- 


Mb,2 

iVsNia[M2(r“)] =A^ J mB)f 
Mb,i 


dM^ 


( 21 ) 


where M 2 (t“) is the mass of the secondary component, 
a normalization constant and the limits of the integral are 

^ M,(T^) should depend on because small stellar 

populations, characterized by a small amount of massive 
stars, could have an IM F truncated at relatively low masses 
llPloeckinger et ^ . 120141 1. We neglect this effect in this paper. 


2.7.3. Stellar winds and the mass return by 
intermediate-mass stars 

The luminosity produced by stellar winds is taken from re¬ 
sults of the Starburst99 software package ()Leitherer et al.L 
119991 : IVazQuez fc Leithereil 1200511 . We used Padova AGB 
stellar tracks at different metallicities (Z* =0.0004, 0.004, 
0.008, 0.02, 0.05) and used this set of models as a library 
to obtain wind luminosities for each SSP, depending on its 
mass and average metallicity. The energy transfer efficiency 
of the wind power is set to 100 %. 
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The mass-return rate (per unit volume) by winds of 
stars in the (1.0 — 8.0) Mq mass range is calculated as 


5: = n.pl [M4r“)], 


(28) 


where the number of stars with mass 1.0 < M* < 8.0 Mq 
dying per unit time is expressed as 


72.* = 




dr 


0 


if t e [r(8.0 Mq) : r(1.0 Mq)], (29) 

if t 9 [t(8.0 Mq) : t(1.0 Mq)], 


and the mass (per unit volume) returned by dying stars p\ 
in the form of a specific element i includes both the true 
nucleosynthesis yi elds of C, O, and N and the pre-existed 


contribution (see iRecchi et al . 20 0 itl. The nucleosynthe¬ 


sis yi elds are taken from van den Hoek &: Groenewegenl 
(1199711 . 


where Zg is the current metallicity of gas from which stars 
form. Equation (1^ satisfies the expected asymptotic be¬ 
haviour of the stellar metallicity, which approaches that of 
the gas (Z*®'^ —>• Zg) during a massive instantaneous burst 
(AMs —t oo), but remains almost unchanged (2'"®"' —Z*) 
in the quiescent phase AMg —0. We note that step 3 from 
Section needs not be applied to the stellar metallicity, 
because this step takes into account the overall aging of 
stellar populations. 


3. Cooling and heating rates 

We assume that the gas is in collisional ionization equi¬ 
librium and calculate the gas cooling rates using the cool¬ 
ing functions presented in Bbhringer fc Henslerl (|l989ll and 
iDalgarno fc McCravl (Il972ll . 


2.7.4. Star-formation rate 

The phase transformation of gas into stars is controlled by 
the following SFR per unit volume 


V = 


f (pg^epw) 

rp'2 

-^n 

* max(Tff,Tc) 

T^+T^ 

lo 



if rig > He 


Her < He 


if 


and V 
or V • > 0 


where rig and T are the gas number density and tempera¬ 
ture, respectively, V • Dg is the gas velocity divergence, and 
^epw is the cold plus warm gas fraction (see Section 13.41 for 
detail). The exponent 1.5 in Equation (1501) is expected if 
the SFR is proportional to the ratio of the local gas volume 
density to the free-fall time, and the free-fall and cooling 
times are calculated as ts = ^37r/(32Gpg) and Tc = e/A, 
respectivel}0. The term in brackets applies the adopted 
temperature dependence by quickly turning off star for¬ 
mation at temperatures greater th an To = 10^ K. For 
consi stency with other studies (e.g. ISoringel fc Hernauisd 
I 2 OO 3 II . the critical gas number density above which star 
formation is allowed is set to 1.0 cm“^ (see discussion in 
Secti on 101). Finally, the n ormalization constant c* is set to 
0.05 (jSt^son et al'l . [200^ . 

2.8. Stellar metallicity 

The heavy element yields by SNell and also by low- and 
intermediate-mass stars are metallicity dependent. We cal¬ 
culate the stellar metallicity Z* using a procedure similar 
to that applied to the mean stellar age in Section 15751 First, 
we solve for the Lagrangian equation for Z* to account for 
bulk motions 


3.1. Low-temperature cooling 

For the low-temperature cooli ng T < 10^ K, we adopted th e 
cooling rates as described in IDalgarno fc McCravl (|1972ll . 
The chemical elements that are taken into account are H, 
O, C, N, Si, and Fe. For high fractional ionization, cooling is 
mostly provided by collisions of neutral and singly-ionized 
^8 ^ ( 3 (jpiis with thermal electrons. We have assumed that most of 
oxygen and nitrogen is in neutral form. The following equa¬ 
tions describe the cooling efficiency Lg{Xi) (in erg cm^ s“^) 
due to collisions of element Xi of number density with 
thermal electrons of number density Ug. 


- low 


L. 


10-18e-413/T 


- low 


(c+) = (^ 7.9 X 

+ 3.0 X io-i7e-6i90o/T^ ^ 

(Si+) = r- 1/2 g X 

+ 3.0 X 10-17e-63600/T^ ^ 

(Fe+) = ( 4.8 X iQ-^»e-2e94/T 

+ 1.1 X 10-18 [e-884/7’ + 1.3e-98i/^ 


(33) 


(34) 


LriO) = 


7.8 X io-i8e-3496/T^ ^ 
1.74 X 10-24^1/2 
(1 - 7 . 6 T- 1 / 2 ) e-228/^ 


(35) 


0.38 


(1 - 7.7T-i/2j e-326/T 

-b 9.4 X io-232’i/2e-227oo/r^ 
^low(N) = g _2 X 10 - 22 t 1/2 

X (1 - 2.7 X 10 -®t 2 ) e- 27 ™/r. 


dZ* _ 5Z* 7 _ n 


(31) 


(36) 


(37) 


Then, we update the stellar metallicity in every grid cell to 
account for the production of metals due to star formation 
using the following equation 


^new _ 


MsZ* -b AMs Zg 
Ms -b AMs 


(32) 


^ Strictly speaking, this is only true when Tc > rg. Otherwise, 
the density exponent should be set to l.O. 


The cooling rate per unit volume A^°'" (erg cm-8 s”!) due 
to collisions of element Xi with thermal electrons is then 
calculated as L^°'"{Xi)ning. Cooling due to collisions of 
hydrogen atoms with thermal e l ectron s, as tabulated in 
Table 2 of IDalgarno fc McCravl (Il972f) is also taken into 
account. 

For low fractional ionization, collisions of C’*', O, Si’*', 
and Fe+ with neutral hydrogen may become important con- 
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tributors to the total cooling budget. The following equa¬ 
tions describe the corresponding cooling efficiency for Si’*' 
and Fe+. 


Lg"'(Sr) 


'(Fe+) 


= 7.4 X 
= 1.1 X 10 


10-23e-413/T 


-22 


j'e-554/T ^ l_4e-961/T^ ^ 


(38) 

(39) 


The cooling rate per unit volume due to collisions of el¬ 
ement Xi with hydrogen atoms of number density tt-h is 
then calculated as mnH. The cooling efficiencies 

due to collisions of O and C~*~ with neutral hydr ogen are 
taken from Table 4 of iDalgarno fc McCravl (|1972f) . 

To summarize, the total cooling rate at the low- 
temperature regime is calculated as 


A 


low 


erg 



Log temperature (K) 


Log temperature (K) 


E 


n^nnCcpw /ion + Ah"'(^*) ] , (40) 


Fig. 2. Cooling efficiency of gas as a function of tempera¬ 
ture for the ionization fraction /ion = 1.0 but varying metal- 
licity (left panel) and for the gas with solar metallicity but 
varying ionization fractions /ion (right column). 


where summation is done over six elements: H, O, C, N, Si, 
and Fe, where applicable. Here, ^cpw is the cold plus warm 
gas fraction, and /ion = ''T-e/uH is the ionization fraction, 
the value of w hich for a given gas temperature Tg is taken 
from table 2 of ISchure et al.l (|2009ll for the collisionally ion¬ 
ized gas in thermal equilibrium for Tg > 10^ ® K. At lower 
temperatures, /ion is set to a constant value of 0.01. 

We note that most of the considered elements (C, Si, 
Fe, and O) have relatively low critical densities Ucrit at 
which spontaneous and collisional de-excitations become 
comparable for some transition levels. Cooling due to these 
elements at tt-h ^ ricrit is proportional to niUcrit rather 
than to riiriH. We take this into account by multiplying 
Equation (HiH) with rirrit / (rir.rit + nn) a nd taking ricrit from 
table 8 in lHollenbach fc McKed (| 198911 . We not that at a 
low numerical resolution the critical densities may never be 
reached. 


In particular, the left panel has /ion = 1 but varying metal¬ 
licity and the right panel has the solar metallicity but vary¬ 
ing /ion. A strong dependence on the metallicity and ion¬ 
ization fraction for /ion ^ 10~® is evident in the figure. 

3.3. Gas temperature 

In order to use the cooling rates described above, one needs 
to know the gas temperature T. The latter is calculated 
form the ideal equation of state P = pgTZTg/p, where TZ is 
the universal gas constant, using the following expression 
for the mean molecular weight 

Pi 4“ ^e/ion^H 

—"7 I (43) 

/ . 4^2 -f /ion^^H 


3.2. High-temperature cooling 

For the solar-metallicity plasma with temperature 
T > 10"^ K, we h ave a dopted the cooling rates of 

iBohringer fc Hensled (Il989tl . which include the line emis¬ 
sion for most abundant elements and continuum emission 
due to free-free, two-photon, and recombination radiation. 
We have considered nine chemical elements (H, He, C, N, 
O, Ne, Mg, Si, and Fe), the total (line plus continuum) 
co oling efficiencies T^‘^^ ( X,;) fo r which are plotted in Fig. 2 
of IBohringer fc Hensled (|1989ll . These cooling efficiencies 
are normalized to neUn and hence need to be multiplied by 
the ionization fraction /ion = u-e/nn to retrieve the values 
that can later be used in our hydro code. The total cooling 
rates in the high-temperature regime is then calculated as 




erg 


Lem"’ sJ 


^^^high(^^)/.^n (nnCcpw)^^ 


(41) 


where is the current number density of element i , is 
its number density at the solar abundance, and summation 
is performed over nine elements: H, He, C, N, O, Ne, Mg, 
Si, and Fe. 

Figure [2] shows the cooling efficiency laid out in the 
following form 


V — IfionL^nX,) + Al?'"(^z)] + (42) 

^^ ry.zj 


where rUe is the electron mass and the summation is per¬ 
formed over all nine chemical elements considered in this 
study. 

3.4. Avoiding the overcooling of SN ejecta 

When the supernova energy is released in the form of the 
thermal energy, there may exists the risk of overestimating 
the gas cooling, which may lead to radiating away most of 
the released thermal energy. This occurs due to the fact 
that the cooling rate A is proportional to the square of gas 
number density, which, for a single-phase medium, may be 
significantly higher than what is expected in the case of a 
multi-phase medium for the hot supernova ejecta. 

To avoid the spurious gas overcooling, we make use of 
the following procedure. We keep track of the hot SN ejecta 
by solving for the following equation for the hot gas density 
Pg,h 


^ + V-ipg,^Vg)=S-^, (44) 

”^c,h 

The second term on the r.h.s. accounts for the phase trans¬ 
formation of the hot gas into the cold plus warm (CpW) 
component, where Tc^h = e/Ah denotes the cooling time 
of the hot component and the cooling rate Ah is cal¬ 
culated assuming the characteristic cooling efficiency of 


7 








































Vorobyov et al.: Stellar hydrodynamics and evolution of dwarf galaxies 


5 X 10“^^ ergs cm^ of the hot component with the 
typical temperature of 10® K (see Fig. [2] and Section |6]0. 

Equation (1441) is solved along with the system of equa¬ 
tions o-® describing the evolution of the gas disk and the 
CpW gas fraction is calculated as .fcpw = (Pg —Pg.h)/Pg- We 
then modify the cooling rates and the SFR by multiplying 
the total gas density pg and the hydrogen number density 
riH in Equations (HiH) and (1411) and also in Equation (1501) 
with ^cpw to retrieve the density of the CpW component. 
This procedure scales down the cooling rates and the SE 
rate according to the mass of the CpW component in each 
computational cell and helps to avoid overcooling and un- 
physically high SF rates. 


3.5. Heating processes 
3.5.1. Cosmic ray heating 


The cosmic ray heating efficiency is the product of the 
cosmic-ray ionization rate ^cr and the energy input per ion¬ 
ization AQrr- We adopt £nr = 4 X 10 ~^^ s“^ based on the 
work of I Van Dishoeck fc Black! (|1986fl . The value of AQcr 
is approximately 20 eV. The resulting heating rate due to 
ionization of hydrogen and helium atoms is 


Fcr = 1.34 X 10"2^Go (0.48nH + O.SuHe) 


1 ergs 1 
Lcm'^ sJ 


(45) 


w here the coefficients 0 . 48 an d 0.5 are taken from table 6 
in iHollenbach fc McKeel (Il989ll . The quantity Gq represents 
the interstellar radiation held (IRF) normalized to the solar 
neighbourhood value, i.e. we make here the usual assump¬ 
tion that the cosmic ray intensity scales with the IRF. We 
calculate its value as Go = SFRioc/SFRmw, where SFRioc 
is the local SFR in our model s and SFRmw = 1-0 yr~ ^ 
is the SFR in the Milky Way (|Robitaille fc Whitnevl . l20i^ . 


3.5.2. Photoelectric heating from small dust grains 

The photoelectric emission from small grains induced by 
FUV photons can substantially contribute to the ther¬ 
mal balance of the interstell ar medium at tempera tures 
^ 10^ K. We use the results of iBakes fc Tielen^ (1 199411 and 
express the photoelectric heating rate per unit volume as 

Fph = 1.0 X lO’^'^CphGonHCd [ergs cm“®s“^] . (46) 

The photoelectric heating efficiency Cph is calculated from 
Figure 13 of Bakes fc Tielens. We have also introduced 
the coefficient = (fcfd/fcfHi)/(fc^d/fcfHi )0 to take into 
account a (possibly) different dust to gas mass ratio 
(Md/Mni) in DGs. The latter is found using the follow¬ 
ing relation between the oxygen abundance (by number 
density) and the dust-to-g as mass ratio in dwarf galaxies 
(|Lisenfeld fc Ferrari 1199^ 

12 + l„g(2)=O.5201„g(^), (47) 

where the normalization constant Gd = —5.0 is derived 
from the corresponding quantities of the Small Magellanic 
Cloud. For the solar dust to gas mass ratio we adopt 
(Md/Mni)© = 0 . 02 , typical for the interstellar medium in 
the Solar neighbourhood. 

In a more accurate approach, Ah can be directly calculated 
using equation m 


4. Solution method 


The equations of gas and stellar hydrodynamics are solved 
in cylindrical coordinates with the assumption of axial 
symmetry using a time-explicit (except for the gas inter¬ 
nal energy update due to gas cooling/heating), operator- 
split approach similar in methodology to the ZEUS code 
(IStone fc NormanL Il992f) . We use the same finite-difference 
scheme but, contrary to the ZEUS code, we advect the in¬ 
ternal energy density e rather than the specific one e/pg 
because the latter may give rise to un desired oscillations in 
some test problems ('e.g Iciark5l201flh . Below, we briefly re¬ 
view the numerical scheme. The test results of our author’s 
implementation of the code are provided in the Appendix. 

The solution procedure consists of three main steps. 
Step 1: we update the gas and stellar components due to 
star formation and feedback. The corresponding ordinary 
differential equations describing the time evolution of pg, pi, 
PgVg, e, ps, PsVs, and Ily due to the source and sink terms 
involving 22, 22®, S, 5®, and F* in equations ®-® and 
® are solved using a first-order explicit finite-difference 
scheme. In Step 1, we also calculate the total gravitational 
potential of the gas and stellar components by solving for 
the Poisson equation m using the alternative d i rectio n 
implicit method described in iBlack fc Bodenheimeil (|1975ll . 
The gravitational potential at the outer boundaries is cal¬ 
culated us ing a multipole expansion formula in spherical co¬ 
ordinates (iJacksonl Il975ll . To save on computational time, 
we invoke the gravitational potential solver only when the 
relative c hange in the total (gas pl us stellar) density exceeds 
10“® fsee lStone fc Normanl 11993. for details). 

Step 2: we solve for the gas hydrodynamics equa¬ 
tions 0 -® excluding the source and sink terms that have 
been taken into account in the previous step. The solu¬ 
tion is split into the transport and source substeps. During 
the former, advection is calculated u sing a third-order ac¬ 
curate piecewise-parabolic scheme of IColella fc WoodwardI 
(11984 ) modihed for cylin drical coordinates as described in 
iBlondin fc LufkiiJ (|1993D . During the latter, we update the 
gas momenta PgVg due to gravitational and pressure forces 
and the gas internal energy e due to com pressional heatini 
using a solution procedure described in IStone fc Norma: 

(EmI- 

At the end of Step 2, a fully implicit backward Euler 
scheme combined with Newton-Raphson iterations is used 
to advance e in time due to volume cooling A and heating 
F rates. In order to monitor accuracy, the total change in 
e in one time step is kept below 10%. If this condition is 
not met, we employ local subcycling in a particular cell by 
reducing the time step in the backward Euler scheme by 
a factor of two (as compared to the global hydrodynamics 
time step) and making two individual substeps in the cell 
where accuracy was violated. The local time step may be 
further divided by a factor of two until the desired accuracy 
is reached. The local subcycling help to greatly accelerate 
the numerical simulations as one need not to repeat the 
solution procedure for every grid cell. 

In Step 2, we also add tensor artificial viscosity to equa- 
ti ons ® and ® to smoot h out strong shocks as described 
Stone fc NormanI (I1992D. As point ed out in iTasker et al.l 


(|20(^) and discussed in Q arkel (201flD. the ZEUS code with 
the Ivon Neumann fc RiAtmvei (|l95nf) definition of the ar¬ 
tificial viscosity performs poorly on multidimesional tests 
such as the Sedov point explosion, yielding non-spherical so- 
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lutions that may overshoot or undershoot the analytic one. 
We demonstrate in Section lA.21 that using the tensor artifi¬ 
cial viscosity and the proper choice of the artificial viscosity 
parameter C 2 = 6 can greatly improve the code perfor¬ 
mance. However, contrary to the Stone & Norman sugges¬ 
tion to discard the off-diagonal terms of the viscous stress 
tensor, we demonstrate that in fact these terms should be 
retained and the most general formulation of the artificial 
viscosity stress tensor should be used 


Q 


iVgV • Ug [Vug 
0 


e(V • Ug)/3] if V • Ug < 0 . ^ 

if V • Ug > 0 , 


where e is the unit tensor, I = C2'max^dx) , and dx is the 
grid resolution in each coordinate direction. For the reasons 
of the code stability, all components of the viscous sterss 
tensor Q should be defined at the same position on the grid, 
i.e., at the zone centers (contrary to the ZEUS code which 
defines the off diagonal components of rank-two tensors at 
the zone corners). The Courant limitations on artificial vis¬ 
cosity may be substantial in regions with strong shocks. 
Therefore, we apply local subcycling (in the same manner 
as it is done with the internal energy update due to heat¬ 
ing/cooling) when the local viscous time step becomes a 
small fraction of the global hydrodynamics time step. 

In step 3, we solve for the stellar hydrodynamics equa¬ 
tions dH)-® excluding the source and sink terms that have 
been taken into account during step 1. The solution proce¬ 
dure is similar to that for the gas hydrodynamics equations, 
since both systems are defined essentially on the same grid 
stencil—a powerful advantage of the stellar hydrodynamics 
approach. All components of the velocity dispersion tensor 
are defined at the cell centers to be consistent with the def¬ 
inition of the stellar artificial viscosity tensor Q. The latter 
is defined by analogy to its gas counterpart with the gas 
density and velocity in equation (H51) substituted with the 
stellar density and velocity. Additional terms — V • Q and 
-Qik ■ (Vus)fcj - Qjk : iVvs)ki have been added to the 
r.h.s. of equations © and respectively, to take into 
account the viscous stress and heating due to artificial vis¬ 
cosity. 

Finally, the global time step for the next loop of simu¬ 
lations is calculated using the Courant condition modified 
to take into account star formation and feedback. In par¬ 
ticular, two time step delimiters of the form 

Ats.f. < (star formation), (49) 

Atf.b. < —(stellar feedback) (50) 

are added to the usual Courant condition to avoid obtain¬ 
ing negative values of the flow variables. The safety coef- 
hcient C 3 is set to 0.5. Some modification to the classic 
Courant condition is required for the stellar hydrodynam¬ 
ics part because the stellar velocity dispersion tensor 11 is 
anisotropic in general. To calculate the stellar hydrodynam¬ 
ics time step, we use the following form 


Afgt = min 


dx 

|us| -I- maxHii 


(51) 


where max H^^ is the maximum diagonal element of the di¬ 
agonalized velocity dispersion tensor and minimum is taken 
over all grid zones. 



Fig. 3. Initial radial profiles of the gas surface density Eg 
(right) and Toomre Q-parameters (left) in models with dif¬ 
ferent rotational support against gravity as indicated by the 
a-parameter. The horizontal dotted lines define the critical 
values of Q-parameter below which star formation is sup¬ 
posed to occur (see the text for more details). 


5. The chemo-dynamical evolution of model dwarf 
galaxies 

Our numerical code has been extensively tested using a 
standard set of test problems as described in the Appendix. 
However, even a perfect performance on the standardized 
tests does not guarantee meaningful results on real prob¬ 
lems when all parts are invoked altogether. 

Therefore, in this section, we consider the evolution of 
several representative models of DCs in order to test the 
overall code performance and to demonstrate the simula¬ 
tion methodology and recipes described in the previous sec¬ 
tions. In particular, we wish to emphasize the effect of stel¬ 
lar dynamics and initial conditions on the simulation out¬ 
comes. Moreover, we will show the results of simulations 
with different degrees of rotational support against gravity 
parameterized by the a-parameter defined as a = v^/vdrc^ 
where and Ucirc are the rotation and circular velocities, 
the latter being the maximum velocity calculated from the 
assumption that all support against gravity comes from ro¬ 
tation (i.e., zero gas pressure gradients). We consider three 
different degrees of rotational support: a =0.25, 0.5, 0.9. 
The resulting initial geometries vary from flat-like for large 
values of a to roundish for small a. 

Our initial configuration consists of a self-gravitating 
gaseous disk submerged into a fixed DM halo described by a 
cored isothermal sphere fe.g. lSilich h Tenorio-TagMl2001ll . 
We focus on models with a DM halo mass of 10® Mq and 
initial gas temperatur es of 10^ K. We apply the solution 
procedure described in IVorobvov et al.l (l2012ll to construct 
self-gravitating equilibrium configurations for different val¬ 
ues of a. The initial gas metallicity Zg is set to 10’ -"^ 0 , 
the ratio of specific heats 7 to 5/3, and the redshift is set to 
zero. The number of grid zones is 1000 in each coordinate 
direction, resulting in the effective numerical resolution of 
8 pc everywhere throughout the computational domain. We 
impose the equatorial symmetry at the midplane and the 
reflection boundary condition at the axis of symmetry. At 
the other boundaries, a free outflow condition is applied. 

Figure |3] presents the initial gas surface density Eg (left 
panel) and the Toomre Q-parameter (right panel) as a func¬ 
tion of the midplane distance in the a = 0.25 model 1 (dash- 
dotted lines), a = 0.5 model 2 (dashed lines), and a = 0.9 
model 3 (solid lines). Evidently, the a = 0.25 model is char¬ 
acterized by the most centrally condensed gas distribution. 
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Table 1 . Model parameters 


Model 

a 

Mg(< 1 kpc) 
Mq 

Mg(< 3 kpc) 
Mq 

Mg (total) 
Mq 

self-consistent 
initial model 

stellar motion 

1 

0.25 

2.1 X 10^ 

2.9 X 10^ 

3.5 X 10' 

yes 

yes 

2 

0.5 

2.8 X 10'^ 

4.5 X 10^ 

6.4 X 10^ 

yes 

yes 

3 

0.9 

1.1 X 10® 

4.1 X 10® 

1.0 X 10® 

yes 

yes 

4 

0.5 

2.8 X 10'^ 

4.5 X 10^ 

6.4 X 10^ 

yes 

no 

5 

0.5 

1.0 X 10'^ 

4.5 X 10^ 

1.2 X 10® 

no 

yes 


The horizontal dotted line marks a critical Q-parameter of 
1.5 below which star formation is supposed to occur ac¬ 
cording to_the stabiliti;_anal^;sis_ofself-gra^ating disks by 
iToomrel (|1964ll and iPolvachenko et all (|1997f) 0. Evidently, 
the a = 0.25 and a = 0.5 models are strongly inclined 
to star formation in the inner 0.5 kpc, while outside this 
radius the conditions are not favourable due to a rather 
steep drop in the gas surface density with distance. On the 
other hand, the a = 0.9 model is characterized by a much 
shallower profile of Eg, resulting in a significantly larger 
galactic volume prone to star formation. The parameters 
of our models are listed in Table [TJ where Mg(r < 1 kpc), 
Mg(r < 3 kpc), and Mg (total) are the gas masses inside a 
spherical radius of 1 kpc, 3 kpc and the total gas mass in 
the computational domain of 8.0 x 8.0 kpc^. 

Evidently, the gas mass increases with increasing a, 
which is explained by the fact that models with higher ro¬ 
tation can support a higher gas mass against gravity of 
the DM halo. We note that the maximum gas mass is also 
limite d by the action of gas self-gravity (|Vorobvov et al.l . 
I2OI2II . As the a = 0.5 model demonstrates, switching off 
the gas self-gravity leads to a factor of two increase in the 
total gas mass. 

When constructing the initial gas density distribution in 
models 1, 2, and 3, we have chosen the same “seed” value for 
gas number density in the galactic center equal to fO cin“^. 
This was done in order to make the initial volume densi¬ 
ties and temperatures in the galactic center similar in each 
model. We could have decreased the seed value in models 2 
and 3, trying to adjust the total mass in these models to 
that obtained in model f, but then the initial conditions for 
star formation would also have changed. In addition, it was 
not clear which mass to take as the reference value, because 
different values are obtained when integrating over different 
radii due to variations in the radial density profiles. 

5.1. The role of a 

In Figure m we show a comparison of the gas volume den¬ 
sity (pg) distribution in three model galaxies differing for 
the value of a. Four representative time snapshots are taken 
for each model. The a=0.9 model 3 is characterized by the 
largest degree of flattening, due to the large centrifugal 
forces. A notable flaring in the initial gas density distri¬ 
bution of model 3 can be seen in the upper-right panel. 
This means that the vertical distribution at i? = 0 is char¬ 
acterized by a very steep density and pressure gradient. 
The over-pressurized superbubble created by the stellar 
feedback preferentially expands along this direction and a 
bipolar outflow is soon created. The external parts of the 

® For a more detailed a nalysis of these star formation criteria 
see IVorobvov et all ll2ni2fl . 


Radial distance (kpc) Radial distance (kpc) Radial distance (kpc) 



0.01 0.1 1 10 0.01 0.1 1 10 0.01 0.1 1 10 



0.01 0.1 1 10 0.01 0.1 1 10 0.01 0.1 1 10 


Radial distance (kpc) Radial distance (kpc) Radial distance (kpc) 

Fig. 5 . Gas surface density as a function of radius for the 
same models and at the same moments in time as in Figure 
SI In the upper panels the initial density distribution is also 
indicated (dashed lines). The dotted lines in the middle 
column also show the distribution for the a = 0.5 model 
without stellar motion (see Section SHI . 

disk are not affected by the bipolar outflow and even af¬ 
ter t = 150 Myr the gas density of the disk at radii larger 
than 1 kpc is quite large. Some of the material in the bipo¬ 
lar outflows falls back onto the disk. As a result, the frac¬ 
tion of pristine gas lost by this model is relatively small, in 
spite of the fact that the galactic outflow is very prominent 
and affects a large fraction of the computational domain. 
However, the freshly produced metals can be easily chan¬ 
nelled out of the galaxy. 

At the other extreme is the a=0.25 model 1. Here, the 
initial gas density configuration is almost spherical and the 
pressure gradient is only slightly steeper in the vertical di¬ 
rection than in the horizontal one. A bipolar outflow soon 
develops but there is significant gas transport also along 
the disk and, after 100 Myr, all the gas in the central part 
of the galaxy (up to a radius of almost 5 kpc) has been 
completely blown away. 

This effect is illustrates in Figure [5l where the gas 
surface density as a function of radius is plotted for our 
three reference models at the same moments in time as 
in Figure 0] After 100 Myr there is almost no gas left in 
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Fig. 4, Gas density distribution of three model galaxies differing for the amount of rotational support: the a=0.25 model 1 
(left column), the q:= 0.5 model 2 (middle column), and the a=0.9 model 3 (right column). The time passed from the 
beginning of numerical simulations is indicated in the leftmost image of each row. The density scale (in cm“^) is shown 
in the upper strip. 


the inner 3 kpc in the a=0.25 model 1, whereas a cen¬ 
tral hole in the q;= 0.9 model 3 is of a much smaller ra¬ 
dius, hardly exceeding 0.2 kpc, and the gas distribution at 
R > 2 kpc is nearly undisturbed. This different galactic 
evolution, depending on the initial distribution of gas, has 
been already described in the literature (see in particular 
iRecchi fc Hensleil l2013^ , but here we show it in a more re¬ 
alistic context of self-gravitating initial gas configurations. 

A direct way to appreciate the role of a on the evolution 
of our model galaxies is by calculating the ejected gas mass 
for various models. Figure [5] presents the fraction of ejected 


gas mass vs. time in models 1-3. We find this value by cal¬ 
culating the amount of gas that leaves the computational 
domain (8x8 kpc^) with velocities greater than the es¬ 
cape velocity at the computational boundary. The fraction 
of ejected mass is very large in the a=0.25 model, because 
of the blow-away described above, but it is very small for 
the q;=0.9 model because only a small fraction of the gas 
originally present in the disk is involved in the galactic out¬ 
flow. The dash-double dotted line presents the fraction of 
ejected mass in the a = 0.5 model 2, in which the motion 
of star was artihcially suppressed (see Sect. 15.21 for details.) 
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Fig. 6. Fraction of ejected gas mass as a function of time for 
the models with different values of a. The dashed-double 
dotted line refers to a model in which the stellar motion 
has been artificially suppressed (see Sect. 



Fig. 7. Star formation rates vs. time in three models dif¬ 
fering in the value of the centrifugal support a. 


The different gas dynamics of the three reference mod¬ 
els shown in Figured has dramatic consequences also for 
the star formation histories. Figure [7] presents the star for¬ 
mation rate vs. time. Evidently, the strong star formation 
feedback in the q;=0.25 model shuts off star formation af¬ 
ter ^ 30 Myr. The a=0.5 model behaves similarly, although 
star formation is shut off slightly later. On the other hand, 
the a=0.9 model is able to sustain star formation for a 
longer time, mainly because the disk still remains gas-rich 
and the conditions for the star formation to occur are ful¬ 
filled in many grid cells. 

Figure m presents the distribution of the stellar volume 
density ps as a function of time and for the same models 
and at the same time instances as in Figure 01 In addi¬ 
tion, the right-hand column presents the volume density of 
high-mass stars p^. We can clearly see that the stellar dis¬ 
tribution in the a = 0.25 model is more homogeneous and 
isotropic than in the a = 0.9 model, wherein many stars 
tend to be located along the polar axis. We note, however, 
that owing to a preferable concentration of gas toward the 
midplane in the a=0.9 model, a significant fraction of stars 
is also distributed near the midplane, constituting a stel¬ 


lar disk. This effect is much less evident for the models 
with q;=0.5 and a=0.25. To quantify this, we calculated 
the fraction of the total stellar mass confined within a ver¬ 
tical height of 1 kpc and 2 kpc. It turns out that in the 
a = 0.25 model these values are 37% and 62%, respec¬ 
tively, whereas in the a = 0.9 model they are 90% and 
94%, indicating a much more compact stellar distribution 
in the higher-a model. Finally, a comparison of the distri¬ 
bution between p^ (the total stellar density) and in the 
q;=0.9 model indicates that stars are preferentially forming 
at low galactic altitudes ^ 1 kpc, but are later spread out 
owing to their own initial dispersion (5 km s“^) and, es¬ 
pecially, owing to significant bulk velocities inherited from 
parental gas clouds. The off-plane star formation is taking 
place predominantly in dense gaseous clumps and on the 
cavity walls. 

5.2. The role of the stellar motion and the initial gas 
distribution 

In this work, we employ the moments of the collisional 
Boltzmann equation to describe the motion of stars. In 
addition, we start our simulations from self-consistent ini¬ 
tial equilibrium conditions that were built taking gas self¬ 
gravity into account (in some studies, the initial configu¬ 
ration is constructed without taking gas self-gravity into 
account, although during the dynamical evolution the gas 
self-gravity may be included). It is instructive to study how 
these two new recipes affect the dynamics of our model 
galaxies. 

In Figure 0] we show the gas volume density in three 
models with q;=0.5 at four representative evolution times. 
In particular, the left-hand-side column shows the evolu¬ 
tion of our reference model 2, already shown in Figure 0] 
The middle column shows model 4, in which the stellar dy¬ 
namics has been artificially suppressed. Stars are allowed 
to form, but they are effectively motionless and are fixed 
to their birthplaces until they die. Evidently, the final gas 
distribution in model 4 is clumpy: the galactic volume is 
filled with moving cometary-like clouds. The increase in the 
number of clumps is caused by a reduced feedback in the 
absence of stellar motion. The moving clumps leave behind 
the newly-born stars so that part of the stellar energy is 
released outside the clumps, which makes it easier for the 
clumps to survive. In the presence of stellar motion, the 
newly-born stars follow the gas clumps (at least for some 
time) due to the fact that they inherit the gas velocity at 
the formation epoch, releasing the stellar energy inside the 
clumps and dispersing most of them. 

The final gas distribution in model 2 is smooth and 
shows no cometary-like structures. At the same time, the 
outflow in model 4 clears a larger hole near the midplane 
than in the reference model. The combination of these two 
effects results in a rather chaotic gas surface density distri¬ 
bution shown in the middle panel of Figure [5] by the dotted 
lines. We also note that model 4 is characterized by a some¬ 
what higher star formation feedback, as is indicated by the 
fraction of ejected gas mass shown by the dot-dot-dashed 
line in Figure 01 

The increased feedback in the model without stellar mo¬ 
tions is likely due to the fact that the energy of both su¬ 
pernova explosions and stellar winds is released in regions 
characterized by a systematically higher temperature and 
lower density than in the reference model. This is because 
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Fig. 8. Similar to Fig.|4j but now for the stellar volume density distribution. The scale bar is in Mq pc The right-hand 
column presents only the massive stars. 


star formation creates regions of slightly lower density com¬ 
pared to the surrounding gas, and stellar winds stir and 
heat up the circumstellar gas. If stars are motionless, then 
their feedback is confined to their birthplaces. On the other 
hand, if stellar motion is allowed, stars can move to neigh¬ 
boring regions, where the density is likely to be higher and 
where stellar winds have not yet heated the gas. The feed¬ 
back will be thus less effecti ve. This conclusion is at vari¬ 
ance with that of lSlvd (|2007fl , who showed that stellar mo¬ 
tion makes the feedback more effective and can even solve 
the overcooling problem. This different outcome is proba¬ 
bly due to the lack of feedback from stellar winds in Slyz’s 
models. 

Figure ITUl presents the stellar volume density distribu¬ 
tions for the same models as in Figure [5] The left-hand col¬ 
umn corresponds to the reference model 2, already shown 
in Figure [5l while the middle column presents model 4, 
in which stellar motions are suppressed. Evidently, models 
with and without stellar motion produce strikingly different 
stellar distributions. The stellar distribution in the refer¬ 
ence model is rather smooth, whereas in the model without 
stellar motions most of the stellar content is confined to 
the midplane and in the inner region. In fact, the fraction 
of the total stellar mass confined within a vertical height of 


0.3 kpc (from the midplane) is 82% for the reference model, 
but this value increases to 98% for the model without stel¬ 
lar motions. The curious stellar stripes visible in the middle 
column are caused by the development of a gaseous super¬ 
shell; many stars tend to form on the walls of this supershell 
and, as a consequence, their position reflects the shell’s ex¬ 
pansion. The lack of stellar motion makes these structures 
permanent. The fact that these structures are not observed 
in real galaxies outlines the importance of a proper treat¬ 
ment of stellar dynamics in numerical simulations of DGs. 

Finally, we investigate the question of how the choice 
of the initial gas density distribution can influence the sub¬ 
sequent evolution of dwarf galaxies. For this purpose, we 
choose a = 0.5 and construct two initial gas density dis¬ 
tributions: one with the gas self-gravity taken into account 
and the other without it. To make the comparison between 
models with and without-self gravity easier, we require that 
the gas mass within the 3-kpc radius be similar in both 
models. The model with the initial equilibrium distribution 
taking gas self-gravity into account is in fact our model 2 
considered above, while its counterpart is further referred 
as model 5 (see Table [1] for model parameters). 

The initial radial profiles of the gas surface density and 
the Toomre Q-parameter in model 5 are shown in Figure |3] 
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Fig. 9. Time evolution of the gas volume density in three models with a=0.5. The l.h.s. column presents the time 
evolution in model 2, already shown in Fig. IT] The middle column shows model 4, in which the motions of stars are 
artificially suppressed. In the r.h.s. column we show model 5, for which the initial equilibrium distribution was constructed 
without taking gas self-gravity into account. Note the different evolutionary times in the latter model. 


by the dash-dot-dotted lines. Not surprisingly, the initial 
gas distribution in model 2 is more centrally condensed than 
in model 5: gas self-gravity creates an additional gravity 
pull towards the galactic centei0, resulting in a denser ini¬ 
tial configuration. For instance, the total gas mass in the in¬ 
ner sphere of 1 kpc in the initial model without self-gravity 
is ~ 10^ Mq, while this value rises by a factor of 2.8 in the 
model with self-gravity. As a consequence, the Toomre Q- 


° The total gas and DM masses are comparable in the inner 
several kpc. 


parameter is initially greater than 2.0 everywhere and star 
formation is effectively suppressed. 

The right-hand column in Figured presents the gas vol¬ 
ume density distribution in model 5. We emphasize that the 
gas self-gravity in this model is neglected only when build¬ 
ing the initial distribution; in the subsequent evolution the 
gas self-gravity is however turned on again. Evidently, the 
evolution of model 5 is notably slower owing to the fact 
that star formation is initially suppressed. The gas cools 
and contracts and only after 120 Myr the conditions for star 
formation become favourable in the galactic center, whereas 
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a vigourous star formation feedback is seen in model 2 al¬ 
ready after 50 Myr of the evolution. The initial slow start 
in model 5 is consistent with a long dynamical timescale 
in this model. For the mean gas volume density in the in¬ 
ner 1 kpc of 4 X 10“^® g cm“^, the corresponding free-fall 
time in model 5 is approximately 110 Myr, which is consis¬ 
tent with the time epoch of the first star formation episode 
in this model. Due to significant thermal support and ra¬ 
dially increasing dynamical timescale (density drops with 
radius), it takes several dynamical timescales to adjust to a 
state similar to that of model 2. By the end of simulations 
in model 5 (800 Myr), the outflow is still not as strong as 
in its counterpart model 2 and the central gap has not yet 
developed. 

The right-hand column in Figure [TUI presents the stellar 
volume density distribution in model 5. A visual inspection 
of the l.h.s. and r.h.s. models reveals that the stellar distri¬ 
butions in models 2 and 5 are in general similar (but note 
the difference in ages 150 vs. 800 Myr), though the stars 
in the latter model are distributed somewhat more irregu¬ 
larly. By the end of numerical simulations, the total mass 
of formed stars in the 800-Myr-old model 5 is still a factor 
of 4 lower than in the 150 Myr-old model 2, highlighting 
a much slower evolution in the model that starts from a 
non-self-consistent initial configuration, i.e. without initial 
gas self-gravity. 

6. Discussion and model caveats 

In this work, we did not aim at a detailed comparison of our 
numerical simulations with observed DGs, which is our task 
for subsequent studies. Instead, we have chosen three typi¬ 
cal models for DGs with different rotational support against 
gravity to test the code performance on global simulations, 
and not just on a predefined set of classic test problems. The 
current version of the code assumes axial symmetry and ne¬ 
glected chemical reactions, which leaves us a wide door for 
future improvements. At the same time, we have managed 
to develop an efficient mathematical recipe for describing 
the phase transition from stars into gas using the continuity 
equation for the mean stellar age (see Section 2.6). 

As already mentioned in the text, some simulation pa¬ 
rameters have been chosen without a detailed parameter 
study or without an in-depth analysis. We avoid such a de¬ 
tailed study for the sake of clarity and in order to make 
the paper more concise and readable. A deeper parameter 
study will be performed in follow-up papers. It is never¬ 
theless important to point out the potential caveats in the 
choice of some parameters. 

— Csw (see Seet. 2.7.3.): Massive and intermediate-mass 
stars expel a significant fraction of their mass by stellar 
winds. Although their lifetimes shorten with mass, the 
mass-loss rates increase disproportionately, so that for 
solar abundances the total wind energy of most mas¬ 
sive stars exceeds that of SNell ov er the s t ellar l ifetime. 
For radiation-driven HII regions, iLasked l|1967f) stated 
that the ISM is heated only by a very low amount 
and derived analytically an efficiency of 1%, because 
the ionizing radiation energy is almost instantaneously 
lost by line emission and only the expansion of HII re¬ 
gions steers up the surrounding ISM. Although the ad¬ 
ditional wind power amounts to almost 10% of the ra¬ 
diation energy, it is expected to deposit a larger frac¬ 


tion as the thermal and kinetic IS M energies. Since in 
most chemo-dynamic al simulations (iTheis et al.L 119921 : 
ISamland et aLl Il997f) the Ly continuum radiation is ap¬ 
plied with an efficiency of 10“^, Cgw is assumed to 100%. 
This value is also used here. 

” esN (see Seet. 12.7.1]) : For the thermalization efficiency 
of the SNII explosion, i.e. the fraction of the explo¬ 
sion energy transferred to gas thermal energy we have 
adopted the value esN = 1-0 for simplicity and in or¬ 
der to avoid further parameter studies. The thermal¬ 
ization efficien cy has been subject to considerable stud¬ 
ies in the past. [Thornton et al.l ([199811 performed single 
ID supernova explosion models and derived esN = 0.1. 
Estimates from cumulative SNII explosions, forming su¬ 
perbubbles, suggest a significant dependence of esN on 
the temperat u re an d density of th e external medium 
(|Recchi et al.l . 120011 : |R ecc H [20T1 . A different ap¬ 
proach, where the con tribution of a whole stellar pop¬ 
ulatio n is considered (jMelioli fc de Gonveia Dal Find 
l2004fl clearly shows that esN is a function of time. We 
have performed a test simulation adopting esN =0.1. 
The extent of the galactic winds was reduced in this 
case, but not by a large factor. It is in fact important to 
remark that with our formulation for the SFR (see Eq. 
[301), a great deal of self-regulation is achieved. In fact, 
the larger esN, the larger the heating of the surrounding 
medium, the higher the probability tha t the star forma¬ 
tion i s quenched for a period of time ('see lKoepuen et al.L 
[19951 for more details on self-regulated star formation). 
Thus, although we plan to make more careful choices of 
esN in follow-up papers, this parameter appears not to 
change drastically the results of the simulations. 

— Tier (see Sect. l2.7.Ul) . The threshold density above which 
star formation is allowed, rici = 1 cm“^ has been 
chosen in compliance with other nume r ical s tudies 
(ISpringel fc Hernauistl . l2003t iRecchi et al.l . I2007H . Also 
in this case, many different values for Ucr have been cho¬ 
sen in the literature and detailed studies on the effect of 
this parameters have been performed. In the framework 
of ACDM cos mological sim ulations. Ucr of 0.1 cm“^ is 
generally used (|Katz[ . [T9^ and its effec t on th e star for¬ 
mation was studied, e.g., bv iKav et al.l ([2002li . but also 
applied to much larger resol utions of smoo t h-par ticle 
hydrodynamics applications ([Stinson et al.l . 120061) . It 
appe ars that an extremely large value of (up to 100 
cm~^ Iciovernato et al.l . l2010li is required to obtain DGs 
with structur al properties simi l ar to the observed ones 
(but see also iPilkington et al.l . l201lli . Quite generally, 
the star formation density threshold should be set based 
on the density at w hich gravitational inst abilities can 
be resolved (see e.g. [Truelov e et al.l . llOOTli . This leads 
for instance IStinson et al.l ([2(113 ) to assume a threshold 
of 9.3 cm“^ based on their adopted resolution. Clearly, 
our choice of Ucr is simplistic and deeper studies are re¬ 
quired in follow-up publications. We notice only that in 
our simulation only a tiny fraction of gas reaches densi¬ 
ties larger than a few tens of cm“^, therefore the adop¬ 
tion of very large values of Ucr would lead to extremely 
localized and sporadic star formation. We notice also 
th at the very large f eedba ck achieved in the simulations 
of iGovernato et aP ([201f)h is certainly able to produce 
modelled DGs with structural properties similar to the 
observed ones, but it is not clear whether the chemical 
properties can be equally well reproduced. A very large 
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Fig. 10. Similar to Figure O but now for the stellar volume density distribution. 


feedback can create galactic winds, able to expel the 
large majority of metals, freshly produced in the galaxy, 
thus at the end of the evolution the model galaxy should 
be almost deprived of heavy elements. 

— Ah (see Sect. 13.41) . The cooling function typical of hot 
gas Ah, which has been used to calculate the cooling 
timescale of the hot SN ejecta, has been estimated based 
on a temperature of the hot ejecta of 10® K. It must 
be admitted that this temperature only holds just af¬ 
ter the SN explosion, but the ejecta cools significantly 
afterwards. Nevertheless, the thermal SN energy has 
been v aried artificially to probe i ts infl uence on galactic 
winds (|Dalla Vecchia fc Schavel l2012ll . We note, how¬ 
ever, that we are in the cooling regime dominated by 


bremsstrahlung, and the temperature dependence of the 
cooling function is not so severe (the cooling is approx¬ 
imately proportional to Vt in this temperature range). 

7. Conclusions 

In this paper, we presented the stellar hydrodynamics ap¬ 
proach for modeling the evolution of DGs, paying special 
attention to the co-evolution of gas and stars. The dis¬ 
tinctive feature of our models is how we treat the stellar 
component. Unlike many previous studies that use vari¬ 
ous versions of the N-body method to compute the dynam¬ 
ics of stars, we employ the stellar hydrodynamics approach 
based on the moments of the collisionless Boltzmann equa- 
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tion, originally introduced by iBurkert fc Hensler (19^), 
and applied an d furth er developed bvlTheis et al l (1992); 
Samland et alJ (Il997ll : Vorobyov fc Theisl ( 20061 2008 ): 


Mitchell et alJ (12(11.11) : iKulikov ( 2014 ). For the first time, 


we provided an effective mathematical recipe for describing 
the death of stars in the framework of stellar hydrodynam¬ 
ics using the concept of a mean stellar age. We performed 
an extensive testing of our numerical scheme using classical 
and newly developed test problems, paying particular at¬ 
tention to the proper recipe for artificial viscosity in order 
to achieve the best agreement with the analytical solution 
for the Sedov test problem. 

We demonstrate the success of our numerical approach 
by simulating the evolution of three DGs differing in 
the amount of rotation described by the spin parameter 
0=0.25, a=0.5, and 0=0.9. In agreement with previous 
studies, models with low values of o are characterized by 
stronger outflows, ejecting a larger amount of gas from the 
galaxy. This is due to the fact that the initial gas distri¬ 
bution in the low-o models is almost spherical, so that the 
gas pushed by the stellar feedback expands almost isotropi¬ 
cally. On the contrary, models with high values of a exhibit 
a disk-like initial gas distribution with a pronounced radial 
flaring. The pressure gradient in the vertical direction in 
these models is very steep, which favours the development 
of bipolar galactic outflows. The transport of gas in the hor¬ 
izontal direction (i.e. along the galactic midplane) is very 
limited, and therefore this model is able to keep most of 
the initial gas bound to the galaxy. The fraction of ejected 
gas (i.e., the fraction of gas that leaves the computational 
domain with a speed higher than the escape speed) is only 
a few percent for the a = 0.9 model, but it raises to 80% 
for the model with a = 0.25. 

At the same time, models with a low value of a pro¬ 
duce a more diffuse stellar population than models with 
higher a. Only models with a > 0.5 reveal the formation 
of a definite stellar disk component. We demonstrate the 
importance of the stellar dynamics by artificially turning 
off stellar motions. In this case, the resulting stellar pop¬ 
ulation is mostly concentrated to the midplane and shows 
stripe-like structures emanating from the galactic center, 
which are not observed in real galaxies. 

Stars start heating the natal gas almost instantaneously 
owing to the effect of stellar winds neglected in many re¬ 
cent hydrodynamical simulations. We found that if stel¬ 
lar dynamics is suppressed, Type II supernovae (the main 
source of mechanical feedback) explode in a medium heated 
and diluted by the stellar winds. In such an environment, 
radiative losses are small and the supernova feedback is ex¬ 
tremely effective. If, on the other hand, stellar motions are 
allowed, they can move from the natal site (where the gas 
has been heated and stirred by the effect of stellar winds) to 
neighboring regions before exploding as supernovae. These 
regions are characterized by lower temperatures and higher 
densities, so that the supernova feedback becomes less effec¬ 
tive. It is thus crucial to correctly simulate stellar dynam¬ 
ics (together with all relevant sources of feedback) in order 
to study the interplay between stars and the interstellar 
medium. 

Unlike many other studies of DGs, we used the self- 
consistent initial gas configuration that was constructed 
takin g the gas self-gravity into account (|Vorobvov et al.l . 
I 2 OI 2 D . This gas configuration is usually more centrally con¬ 
densed and prone to star formation than that constructed 


neglecting gas self-gravity. We compared the evolution of 
our model galaxy starting from these different initial con¬ 
figurations and found that models that start from non-self- 
gravitating initial configurations evolve much slower than 
those that start from self-gravitating ones. In the former, 
a major episode of star formation and the development of 
prominent outflows may by delayed by 0.5-0.8 Gyr as com¬ 
pared to the latter. We note that this is exclusively the 
effect of initial conditions (or, better to say, the effect of 
using non-self-consistent initial gas configurations), because 
the gas self-gravity is turned on in both models once the 
evolution has started. 

In the future, a detailed comparison of our stellar hy¬ 
drodynamics simulations with a more conventional N-body 
treatment of stellar dynamics is desirable to assess the weak 
and strong sides of the Boltzmann moment equation ap¬ 
proach. This project was partly supported by the Russian 
Ministry of Education and Science Grant 3.96I.2014/K and 
a Lise Meitner Fellowship, project number M I255-NI6. 
The simulations were performed on the Vienna Scientific 
Cluster (VSC-2). This publication is supported by the 
Austrian Science Fund (FWF). We are thankful to the 
anonymous referee for a very insightful review which helped 
to improve the manuscript. 
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Appendix A: Testing gas hydrodynamics equations 

Our numerical scheme for solving the equations of gas hy¬ 
drodynamics 0-0 can be tested using a conventional set 
of test problems suitable for cylindrical coordinates. Below, 
we provide the results of six essential test problems, each 
designed to test the code performance at various physical 
circumstances. The code has also performed well on two ad¬ 
ditional tests: the two interacting blast waves and the Noh 


Distance Distance 



Fig. A.l. Sod shock tube problem for the gas density 
(upper-left), velocity (upper-right), pressure (lower-left), 
and specific internal energy (lower-right). The numerical 
solution is shown by open circles, while the analytical one 
is plotted by solid lines. 


(or implosion) problems, but we do not provide the results 
for the sake of conciseness. 

A.l. Sod shock-tube test 

This test is often used to assess the ability of a numerical 
algorithm to accurately track the position of relatively weak 
shock waves and contact discontinuities. Initial conditions 
involve two discontinuous states along the z-axis, with a 
hot dense gas on the left and cold rarified gas on the right. 
More specifically, we set the pressure and density at z G 
[0 — 0.5] to 1.0, while at z G [0.5 — 1.0] the pressure is 0.1 
and density is 0.125. The velocity of a 7 = 1.4 gas is initially 
zero everywhere. 

Filled circles in Figure ET] show the results of the Sod 
shock-tube test computed with a resolution of 200 grid 
zones and C 2 = 4. (the number of zones over which the shock 
is spread), while the solid lines plot the analytic solution at 
t = 0.235. It is evident that the code tracks well the posi¬ 
tions of shock waves and contact discontinuities, which are 
resolved by 3-4 zones, in accordance with the chosen value 
of C 2 ■ Decreasing C 2 results in sharper shocks and contact 
discontinuities but may increase the magnitude of anoma¬ 
lous spikes in the specific energy (bottom-right panel), and 
hence not recommended. Increasing C 2 to 6 produces re¬ 
sults similar to those for C 2 = 4, except that the shocks 
and contact discontinuities are now resolved by 5-6 zones. 

A.2. Sedov point explosion 

The Sedov explosion problem tests the code’s ability to deal 
with strong shocks. In contrast to the Sod shock-tube test, 
in which shock waves were of a pure planar geometry, this 
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Fig. A.2. Sedov point explosion test for the explosion en¬ 
ergy 10®^ erg and background density n = 0.01 cm“^. 
The analytic solution is shown by the red lines, while the 
numerical solutions along the z-axis, r-axis, and diagonal 
z = r are shown by the black, pink, and blue lines, re¬ 
spectively. The top-left panel presents the solution for the 
tensor artificial viscosity with C 2 = 6, the top-right panel— 
for the scalar artificial viscosity and (72=3, the bottom-left 
panel—for the tensor artificial viscosity and C 2 = 3, and 
the bottom-right panel—for the tensor artificial viscosity 
with zero off-diagonal terms and (72=6. 


test involves a spherically symmetric strong shock wave. 
Hence, the point explosion is a powerful test on the code’s 
ability to render spherically symmetric structures on an 
essentially non-spherical grid stencil. 

To initialize the explosion, we consider a cold (T = 
10 K) homogeneous medium with number density n = 
10 “^ cm“^ and inject 10®^ ergs of thermal energy (equiv¬ 
alent to one supernova) into the innermost grid cell. The 
adopted resolution is 300 x 300 grid cells and the size of 
each cell is 1.0 pc in each coordinate direction {z and r). 
We emphasize that we do not inject energy into a central 
sphere comprising a few cells near the coordinate origin—a 
usual practice to alleviate the problem of a non-spherical 
geometry. Instead, we consider the most difficult situation 
when energy is injected in just one, essentially cylindrical 
innermost cell. 

Figure [X)2l compares our derived density distributions 
along the z-axis (black line), the r-axis (pink line), and the 
diagonal z = r (blue line with circles ) with the analytic 
solution given by the red line (|Sedovl . 11959(1 at t=l Myr 
after the explosion. In particular, numerical simulations in 
the top-left panel are obtained using the tensor artificial 
viscosity (AV) as defined in Section (U while the density 
distribution in the top-right panel i s obtained using the 
” class ic” scalar AV as described by iBichtmver fc MortonI 
(|1995r) . It is evident that the tensor AV yields a much better 
agreement with the analytic solution. The expanding shell 
has almost a perfect spherical shape (all three lines showing 
the numerical solution virtually coincide), notwithstanding 
the fact that the energy was injected into a cylindrically 
shaped innermost cell. On the other hand, the Richtmyer & 
Morton scalar AV yields a notable deviation from spheric¬ 
ity along the r and 2 axes, as seen in the top-right panel of 


Figure [221 The Richtmyer & Morton solution also slightly 
overshoots the analytical one. Regardless of the AV pre¬ 
scription used, the shock front is resolved by 2-3 grid zones. 

We want to emphasize that all components of the ar¬ 
tificial viscosity tensor Q should be used, including the 
off-d iagonal terms. Neglectin g the latter, as is done in, 
e.g., (Stone fc NormanI ((l 992(1 . leads to a significant dete¬ 
rioration of the numerical solution. The bottom-right panel 
in Figure E21 demonstrates the effect of zeroing the off- 
diagonal terms in Q. As one can see, the numerical solution 
is skewed, notably undershooting the position of the shock 
front in the r-direction and along the diagonal direction but 
reproducing the shock front rather well in the z-direction. 

A caution should be paid when choosing the AV softener 
(72. As one can notice in Figure (A.21 the test was done 
with different values of C 2 for the tensor and scalar AV, 
6 and 3, respectively. If we choose C 2 = 3 for the tensor 
AV, the position of the shock front notably undershoots the 
analytic solution, as is illustrated in the bottom-left panel of 
Fig. (A.21 Fortunately, once the proper value of C 2 is found, 
the tensor AV algorithm performs well regardless of the 
energy input, background dens ity, o r numerical resolution. 
We demonstrate this in Figure 123] showing the tensor AV 
performance for various choices of the input parameters as 
indicated in each panel and (72 = 6. It is evident that the 
test results are virtually insensitive to a particular choice 
of initial conditions. 

The sensitivity of the Sedov test to a particular choice 
of the AV softener C 2 is a well known probl em of numeri¬ 
cal co des that evolve the internal energy (e.g. (Tasker et al.L 
(2008(1 . Switchi ng to the tota l energy helps to eliminate this 
problem (e.g. (Clar kilMi. However, in numerical codes 
that make use of a staggered grid as our own, with vec¬ 
tor and scalar variables defined at different positions on 
the grid, solving for the total energy equation is not a 
good option as it inevitably requires interpolating between 
the internal a nd kinetic energies to obtain the total one 
da arkd . (2010(1 or applying co rrections after each tim e step 
to conserve the total energy ([Vshivkov et al.( . [2011(1 . In ei¬ 
ther case, this practice, according to our experience, of¬ 
ten comes at a cost of degraded performance on other test 
problems and is not recommended. Instead, we show that a 
proper choice of the artificial viscosity softener can help to 
resolve this problem without resorting to the total energy 
equation. 


A.3. Collapse of a pressure-free sphere 

The gravitational collapse of a pressure-free sphere is used 
to assess the code’s ability to accurately treat converging 
spherical flows in cylindrical coordinates. This test is also 
useful for estimating the performance of the gravitational 
potential solver on dynamical problems (our solver does 
well on the static configurations such as spheres and ellip¬ 
soids). To run this test, we set a cold homogeneous sphere of 
unit radius and density (for convenience, the gravitational 
constant is also set to unity) and let it collapse under its 
own gravity. The analytic solution to this problem exists 
only in the limit of an infini te cloud radiu s, describing the 
collapse of every mass shell ([Hunted . (l962(l . In the cylindri¬ 
cal geometry, however, we consider a cloud of finite radius 
with a sharp outer boundary to preserve the cloud spheric¬ 
ity. As a result, a rarefaction wave develops after the onset 
of the collapse, propagating towards the coordinate origin 
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Fig. A.3. Sedov point explosion test for various explosion 
energies, background densities and numerical resolutions as 
indicated in each panel. The tensor artificial viscosity with 
(72 = 6 is used for every plot. The line meaning is the same 
as in Figure IA.2I 


and nece ssitating complicated corrections to the analytic 
solution (|Trnelove et al.l . [199811 . 

The left panel in Figure Ed] compares the results of our 
numer ical simulation s with the ” uncorrected” analytic solu¬ 
tion of iHunteJ (Il962fl (red line) at (=0.535 (or 0.985 that of 
the free-fall time) when the initial density has increased by 
nearly three orders of magnitude. As in the previous test, we 
plot the numerically derived density along the z-axis (black 
line), r-axis (green line), and the z = r diagonal (blue line 
with circles). The numerical resolution is 200 x 200 grid 
zones. It is seen that the cloud’s outer boundary is some¬ 
what smeared out due to the action of the rarefaction wave. 
The peak density is also about 1% smaller than that pre¬ 
dicted from the analytic solution. On the other hand, the 
code performs well on preserving sphericity of the collapsing 
cloud. In the right panel, we turn off the volume averaging 
corrections applied to the advec tion of the r-momentum in 
the r-direction (as described in iBlondin fc Luf^ l|l993ri l. 
which results in the development of an artificial density 
spike near the z-axis and also in the distortion of a spheri¬ 
cally symmetric collapse. This test demonstrates the utility 
of the volume aver a ging c orrections introduced originally by 
iBlondin fc Lufkin! (1 199311 in the PPA advection scheme in 
non-Cartesian coordinates. 


A.4. Collapse of a rotating sphere 

This test is invoked to assess the code’s ability to con¬ 
serve angular momentum. The setup consists of a isother¬ 
mal (T=10 K) homogeneous sphere of unit radius and den¬ 
sity, which rotates at a constant angular velocity CIq. The 
latter is found from the requirement that the ratio /3 of 
the rotational energy i?rot = 7I2 q/ 2 to gravitational energy 
= 3GM^/{5R) is equal to 5%. Here, I = 2 Moi?V 5 
is the moment of inertia of a homogeneous sphere of radius 
R and mass Mq. The adopted value of /3 = 0.05 is close to 
an upper limit found in rotating pre-stellar molecular cloud 
cores and dark matter halos. 
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Fig. A. 4. Collapse of a pressure-free sphere showing the 
gas density at 0.985 that of the free-fall time. The left/right 
panels show the test results with /without the volume aver - 
aging corrections as described in iBlondin fc Lufl^ (|1993ll . 
The line meaning is the same as in Figure lA. 2 1 



Specific angular momentum, K 

Fig. A. 5. Specific angular momentum spectrum M(K), cal¬ 
culated using equation for a collapsing sphere at 

1.03(ff. The analytic spectrum is shown by the solid line, 
while the test results are plotted by open circles. 


In the absence of any mechanisms for angular momen¬ 
tum redistribution, the mass with specific angular momen¬ 
tum less than or equal to K = ojr^ , where oj is the angular 
velocity at a d istance r from the z-a xis, should be conserved 
and equal to (|Norman et al.l . Il98flll 


M{K) = Mq 


1 - 1 - 


K 


\ 3/2' 


(A.l) 


A deviation from equation (lA.ll) would manifest the non¬ 
conservation of angular momentum in the numerical algo¬ 
rithm. 

The solid line in Figure ED shows the initial mass spec¬ 
trum M{K), whereas the open circles present the mass 
spectrum obtained at time 1.03 x (g = 0.551, where tg is the 
free-fall time for a non-rotating sphere of the same density 
and radius. The numerical resolution is 200 x 200 grid zones. 
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By t = 0.551, the density near the center of the sphere has 
increased by more than three orders of magnitude, but the 
deviation of the mass spectrum from the initial configu¬ 
ration is negligible everywhere except at lowest values of 
K. This deviation is a manifestation of angular momentum 
non-conservation near the rotational axis caused by im¬ 
perfections in the interpolation procedure across the inner 
r—boundary. Nevertheless, the total mass of gas that suf¬ 
fers form the angular momentum non-conservation is just 
below 1% and therefore should not affect notably the global 
evolution. 

A.5. Kelvin-Helmholtz Instability 

The Kelvin-Helmholtz instability (KHI) occurs at the in¬ 
terface between two fluid moving with different velocities. 
In the context of this paper, this instability operates at the 
interface between the hot supernova ejecta and the cold 
swept-up shell an d may cause an accelerat ed disintegration 
of the latter fe.sr. IVorobvov fc Basd . [200^ . Therefore, it is 
essential that the code be able to demonstrate the develop¬ 
ment of the KHI as expected from the linear perturbation 
theory. 

The initial setup involves two fluids moving in oppo¬ 
site directions along the z-axis with the relative velocity 
fz,rei = 0.2cs, where Cs = 0.06 is the sound speed. The inter¬ 
face between the fluids is located at r = 0.5 and the densi¬ 
ties of the inner and outer fluids are p-m = 1.0 and pout=0.1, 
respectively. We use periodic boundary conditions and the 
resolution is 200 x 200 grid zones on the 1.0 x 1.0 compu¬ 
tational domain. 

To trigger the instability, we impose a sinusoidal per¬ 
turbation on the radial velocity fr of the form 

= i5fr sin(27rz/A), (A-2) 

where the amplitude and wavelength of the perturbation 
are Sv,; = Vz^rei/lOO and A = 1/6, respectively. The per¬ 
turbation is applied to five neighboring grid zones on both 
sides of the fluid interface. 

Figure IA.6I demonstrates the growth of the KHI with 
time. There are six growing vortices, in agreement with the 
wavelength of the perturbation A = 1/6, and the growth 
time is approximately 50-60, again in agreement with the 
characteristic growth time of the KHI for our choice of the 
fluid density and velocity 

tkh = = 48.3. (A.3) 

^z,relY PinPout 

We should note here that, for this particular test, gravity 
can be neglected. Therefore, all modes of perturbation, irre¬ 
spective oKtheimwayelenghtsJor wavenumbers) are unsta¬ 
ble (see IVieser &: Hensleii l2007l for stability criteria when 
gravity is taken into account). 

A.6. Overstability of a radiative shock 

A useful test problem for hydrodynamic codes with opti¬ 
cally t hin radiative cooling is described in IStone fc NormanI 
(|l993fl . The test involves setting a unform flow of gas with 
velocity vq against a reflecting wall. An adiabatic reflection 
shocks forms immediately at the wall and propagates up¬ 
stream with velocity Vsh = r’o/3. After approximately one 
cooling time, the shock looses its pressure support and is 



Fig. A. 6 . Development of the KHI at different moments in 
time (see text for more detail). 



Fig. A.7. Shock position versus time of a one-dimensional 
radiative shock during overstable oscillations. The shock 
positions for the flow velocities of 120 km s“^, 130 km s“^ 
and 150 km s“^ are plotted with solid, dash-dotted, and 
dotted lines, respectively. 


then advected back to the wall at Vsh ^ —vq. After reflect¬ 
ing off the wall, the shock again becomes adiabatic and 
propagates outward to repeat the cycle. 

Figure Et] presents the test results for the cooling func¬ 
tion with solar metallicity and ionization fraction /ion = 
10“^. Gas heating is turned off for this test. The initial 
setup is identical to that described in IStone fc Normaiil 
(Il993ll . i.e., the initial gas number density and temperature 
are set to 10 cm“^ and 10^ K, respectively. The position 
of the shock front for three flow velocities of 120 km s“^, 
130 km s“^ and 150 km s~^ are plotted with the solid, 
dash-dotted, and dashed lines, respectively. The shock po¬ 
sition demonstrates regular pulses, as expected. The ampli- 
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tude and period of the observed oscillations are somewha t 
different from those presented in IStone fc NormanI (llQQilll , 
which simply reflect the difference in the adopted cooling 
functions. The upstream velocity of the adiabatic shock is 
close to the expected value, no/3. For instance, in the case 
of no = 150 km s“^, the maximum upstream velocity of 
the shock is 47 km s“^. This test demonstrates the reliabil¬ 
ity of our solution scheme for the gas cooling and heating 
described in Section |31 

Appendix B: Testing stellar hydrodynamics 
equations 

Testing the stellar hydrodynamics equations 
presents a certain challenge since stellar systems are 
described by a stellar dispersion tensor, which may 
be anisotropic, rather than by an isotropic pressure. 
Fortunately, some of the test problems considered above 
can be adapted to stellar hydrodynamics as well. The mat¬ 
ter is that the equations of stellar hydrodynamics can be 
formally made identical to those of gas hydrodynamics for 
a specific set of initial conditions (note that the continuity 
equations are always identical). Indeed, for a one dimen¬ 
sional flow along the z-direction (with zero gravity and no 
star formation), equations d?]) and ([5]) become 

d d d 

^ {PsVs,z) + ^ {PsVs^z • Vs,z) + ^ {Ps<^zz) — (^. 1 ) 

^ (Pscrlz) + Yz 

These equations become identical to the corresponding gas 
dynamics equations for pgV^ and e, if we set PsO'L = P , 
e = (7 — l)ps0-^2, and 7 = 3 (see iMitchell et al.[ 120131 
for details). This fact allows us to directly compare the 
performance of both stellar and gas hydrodynamics code on 
some test problems considered in Section]^ or use analytic 
solutions available from the gas dynamics. 

B.l. Sod shock-tube test 

The initial setup for this test is identical to that considered 
in Secion fOl We run the test along the z-axis and set 
and ps aX z € [0 — 0.5] to 1.0, while at z € [0.5 — 1.0] the 
z-component of the stress tensor is 0.1 and stellar density 
is 0.125. 

Our numerical results for C 2 s = 4 and t = 0.14 (open 
circles) are compared against the analytical solution (solid 
lines) in Figure IB.II It is seen that our numerical scheme 
correctly reproduces the position of shock waves and con¬ 
tact discontinuities in the stellar fluid. On the other hand, 
small-scale oscillations are evident at shock positions, per¬ 
haps reflecting the lack of energy dissipation due to the 
collisionless nature of stellar hydrodynamics. An increase 
in the coefficient of artificial viscosity helps to reduce the 
amplitude of spurious oscillations but simultaneously in¬ 
creases the shock smearing and hence is recommended only 
if these oscillations lead to numerical instabilities. 

B.2. Collapse of a cold sphere and angular momentum 
conservation 

The collapse of a cold stellar sphere, both with and with¬ 
out rotation, have produces very similar results to those of 


Distance Distance 



Distance Distance 

Fig. B.l. Sod shock tube problem for the stellar den¬ 
sity (upper-left), velocity (upper-right), z-component of the 
stress tensor (lower-left), and square of the z-velocity dis¬ 
persion (lower-right). The numerical solution is shown by 
open circles, while the analytical one is plotted by solid 
lines. 


presented in Sections fA.3l and [4.41 and will not be repeated 
here. This ensures that our stellar hydrodynamics part is 
expected to perform well on problems involving gravita¬ 
tional contraction/expansion with and without rotation. 


Appendix C: Testing the chemical evolution 

routines 

Nowadays it is common to find hydrodynamical codes cou¬ 
pled with passive scalars (or other numerical techniques), 
able to follow the evolution of the metals in a simulation. 
However, usually no test is performed to check the numer¬ 
ical accuracy of these routines. Actually, since a very vast 
literature on analytical solutions of the chemical evolution 
of galax i es an d other astr o nomical objects exists (see e.g . 
iTinslevl . 119801: iMatteuc^ . 120011: iHensler &: Recchi 1201OL 
and references therein) it is indeed not difficult to bench¬ 
mark the chemical evolution routines of a chemodynamical 
code. 

The “closed box” model of chemical evolution is based 
on the following assumptions: (i) the system is uniform and 
closed (there are no inflows or outflows), (ii) the IMF does 
not change with time, {in) the gas is well mixed at any time 
and {iv) stars more massive than a certain threshold mass 
Mthr die instantaneously whereas stars less massive than 
Mthr live forever. This last assumption, although very re¬ 
strictive, is fulfilled if we look at specific chemical elements, 
in particular O and the other a-elements, since they are 
mostly synthesized by Type II SNe and the (long-living) 
low- and intermediate-mass stars contribute negligibly to 
their production. 
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In a system made up of only gas and stars (DM does 
not play any role here) we define the quantities 

pCO 

R= [M - Mrem{M)](l){M)dM, (C.l) 

JMthr 

yo = [ Mpo{M)cj){M)dM, (C.2) 

1 - -K JMthr 

where Mtot = Mgas + Mstars, Mrera{M) Is the remnant 
mass after a star of mass M has died and Mpo{M) is 
the oxygen mass newly synthesized by the star of mass M 
{po{M) is the true yield of O). Under the initial condition 
MgasfiO) = Const. = Mtot, MstarsiO) = 0 , Xo(0) = Xo it 
is not difficult to find that 

Xo = ^ = Xo + yo\n{p-^), (C.3) 

where p = Mgas/Mtot- This is the famous analytical solu¬ 
tion of the closed box model, which depends only on the 
gas fraction p. 

We have integrated numerically eqs. IC.ll and 1C.21 by 
usi ng the values of Mrom/M ) and Mpo{M) tabulated 
by IWooslev &: Weav^ (|1995ll . using Z=0.01 Zq as initial 
metallicity. We do not use Z=0 as initial metallicity for 
two reasons: firstly, for Z(0)=0 the cooling is very low 
and the resulting star forma tion is too weak; secondly, the 
yields po{M) tabulated by IWooslev fc Weav^ (| 199511 at 
this metallicity oscillate considerably and that makes the 
calculation of yo by means of Eq. 1C.21 less accurate. By 
properly choosing a reference volume in the computational 
box (namely a volume encompassing all the metals pro¬ 
duced and advected within a time span of 100 Myr) we 
can thus calculate the variation of p with time and, from 
it, the analytical evolution of the oxygen mass fraction Xq 
with time. This is plotted in Fig. IC.ll (solid line), together 
with the direct numerical calculation of Mq/M gas, aver¬ 
aged over the same volume (dashed line). The two curves 
converge only after ^ 30 Myr. This is expected because 
we have to ensure that Type II SNe of all initial masses 
contribute to the production of O and ^ 30 Myr is approx¬ 
imately the lifetime of the less massive SNell. After a time 
of ~ 40 Myr the two curves almost overlap, demonstrating 
that our chemical routines are able to accurately follow the 
overall chemical evolution of a galaxy. 



Fig. C.l. The evolution with time of the mass fraction 
of oxygen as derived from the analytical expression Eq. 
1C.31 (solid line) and from our numerical simulation (dashed 
line). See text for more detail. 
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